Cosmic microwave background with Brans-Dicke Gravity: I. Covariant Formulation 
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In the covariant cosmological perturbation theory, a 1+3 decomposition ensures that all variables 
in the frame-independent equations are covariant, gauge-invariant and have clear physical inter- 
pretations. We develop this formalism in the case of Brans-Dicke gravity, and apply this method 
to the calculation of cosmic microwave background (CMB) anisotropy and large scale structures 
(LSS). We modify the publicly available Boltzmann code CAMB to calculate numerically the evo- 
lution of the background and adiabatic perturbations, and obtain the temperature and polarization 
spectra of the Brans-Dicke theory for both scalar and tensor mode, the tensor mode result for the 
Brans-Dicke gravity are obtained numerically for the first time. We first present our theoretical for- 
malism in detail, then explicitly describe the techniques used in modifying the CAMB code. These 
techniques are also very useful to other gravity models. At last, we compare the CMB and LSS 
spectra in Brans-Dicke theory with those in the standard general relativity theory. Constraints on 
Brans-Dicke model with current observational data is presented in a companion paper [l)] (paper II). 



I. INTRODUCTION 

The Jordan-Fierz-Brans-Dicke theory 0, 0, 0, S @] 
(hereafter the Brans-Dicke theory for simplicity) is a nat- 
ural alternative and a simple generalization of Einstein's 
general relativity theory, it is also the simplest example 
of a scalar-tensor theory of gravity @, H, S USUI- In the 
Brans-Dicke theory, the purely metric coupling of matter 
with gravity is preserved, thus ensuring the universality 
of free fall (equivalence principle) and the constancy of 
all non-gravitational constants. From early on, testing 
the Brans-Dicke theory with CMB anisotropy has been 
considered fl3 ]. However, the usual approach is to use 
a metric-based and gauge-dependent method, i.e. mak- 
ing the calculation with a particular gauge, see. e.g., 

Ref. 0B Hill, [HI- 

The covariant approach to general relativity is an ele- 
gant solution to the "gauge problem" , which has plagued 
the study of linear perturbation in gauge-dependent 
methods since the pioneering work of Lifshitz[l8[. Before 
this problem was recognized, contradictory predictions 
of the behavior of perturbation of Friedmann-Lemaitrc- 
Robertson- Walker (FLRW) cosmologies were made. In 
1980, Bardeen reformulated the metric approach using 
gauge-invariant variables 0] (see also Ref. 0] for a re- 
view on the variables which has been widely used in re- 
cent perturbation calculations). However, as pointed out 
by Ellis [2l| , although the Bardeen variables are related 
to the density perturbations, they are not those perturba- 
tions themselves, since they include metric tensor Fourier 
components and other quantities in cunning combina- 
tions. The physical meaning of Bardeen's gauge- invariant 
variables are not always transparent. As emphasized by 
Hawking [22| , the metric tensor can not be measured di- 



* Electronic address: 
t Electronic address: 



wufq@bao.ac.cn 



xuelei@cosmology.bao.ac.cn 



rectly, so it is not surprising that the variables used in the 
metric-based method do not always have a clear physical 
interpretation. 

The covariant approach to general relativity and 
cosmology has its origins in the work of Heckmann, 
Schucking, Raychaudhuri and Hawking 0, [H, 0] ■ In 
1989, Ellis and Bruni proposed using the spatial gradient 
of matter density (D a p) as the basic variable to describe 
density perturbations [2]| . Subsequently, the cosmologi- 
cal applications have been dev elop ed extensively by Ellis 
and others in recentyears j2pjM H3, H, [U, 0, IH 

EI Sa 

41]. The method also has 

been applied to problems in CMB physics nasi, S3, El. 

Instead of using the components of metric as basic vari- 
ables, the covariant formalism performs a 1+3 split of the 
Bianchi and Ricci identities, using the kinematic quan- 
tities, energy-momentum tensors of the fluid(s) and the 
gravito-electromagnetic parts of the Weyl tensor to study 
how perturbations evolve. The most notable advantage 
of this method is that the covariant variables have trans- 
parent physical definitions, which ensures that predic- 
tions are always straightforward to interpret physically. 
Other advantages include the unified treatment of scalar, 
vector and tensor modes, a systematic linearization pro- 
cedure which can be extended to consider higher-order 
effects (this means the covariant variables are exactly 
gauge- invariant, independent of any perturbative expan- 
sion), and the ability to linearize about a variety of back- 
ground models, e.g. either the FLRW or the Bianchi 
models 0,113. 

A pioneering work in applying the covariant approach 
to Brans-Dicke theory is Ref. [48{ . in which a conformal 
transformation was performed, and calculation was done 
in the Einstein frame. More recently, Ref. 0, 0] chose 
the effective fluid frame, implying D a <p = and tu a b = 0, 
i.e. their foliation selects vorticity-free spacelike hyper- 
surfaces in which <f> = const, hence greatly simplifies the 
calculations. 

The aim of this paper is to construct a full set of co- 
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variant and gauge-invariant linearized equations to cal- 
culate angular power spectra of CMB temperature and 
polarization anisotropies in the cold dark matter frame, 
and to show that the covariant method will lead to a 
clear, mathematically well-defined description of the evo- 
lution of density perturbations. In a companion paper [l| 
(heretofore denoted paper II), we shall apply the formal- 
ism developed in this paper to the latest CMB and large 
scale structure data to obtain contraint on the Brans- 
Dicke parameter. 

In §2, we briefly review the Brans-Dicke theory and 
its background cosmological evolution. The formalism of 
covariant perturbation theory is presented in §3, and the 
numerical implementation and the results are discussed 
in §4. Finally, we summarize and concludes in §5. 

Throughout this paper we adopt the metric signature 

(H ). Our conventions for the Riemann and Ricci 

tensor are fixed by [V a ,V;,]it c = R a b c dU d , where V a de- 
notes the usual covariant derivative, and R a b = R a cb C '■ 
We use d a to represent ordinary derivative. The spa- 
tially projected part of the covariant derivative is de- 
noted by D a . The index notation Ai denotes the in- 
dex string a\...ai. Round brackets around indices denote 
symmetrization on the indices enclosed, square brackets 
denote anti-symmetrization, and angled brackets denote 
the projected symmetric and tracefree (PSTF) part. We 
adopt k — 8ttG and use units with h = c = ks = 1 
throughout. In the numerical work we use Mpc as unit 
for distance. 



II. THE BRANS-DICKE THEORY AND 
BACKGROUND COSMOLOGY 



The Brans-Dicke theory is a prototype of the scalar- 
tensor theory of gravity. One of its original motivations is 
to realize Mach's principle of inertia [5, 6]. It introduced 
a new degree of freedom of the gravitational interaction 
in the form of a scalar field non-minimally coupled to the 
geometry. The action for the Brans-Dicke theory in the 
usual (Jordan) frame is 



eralized to 



1 

16^ 



-.9 



-<t>R+-g» u Vu<l>V u <t> 



s (m) 



(1) 

where <p is the Brans-Dicke field, w is a dimensionless pa- 
rameter, and (S^™) is the action for the ordinary matter 
fields = J d 4 Xy/^g£( m \ Matter is not directly cou- 
pled to (f>, in the sense that the Lagrangian density £' m ) 
does not depend on </>. For convenience, we also define a 
dimensionless field 



Guv 



8ttG 



(3) 



where Tjht' is the stress tensor for all other matter ex- 
cept for the Brans-Dicke field, and it satisfies the energy- 
momentum conservation equation, V^T^"^ = 0. The 
equation of motion for p is 



2lo + 3 



(4) 



here T^ 11 ' 1 — T^ m ^ M is the trace of the energy-momentum 
tensor. The action ([1]) and the field equation |3]) suggests 
that the Brans-Dicke field plays the role of the inverse 
of the gravitational coupling, G e ff(ip) = \ = — , which 
becomes a function of the spacctimc point . 

For background cosmology, we treat the ordinary mat- 
ter as the perfect fluid with the energy density p and 
pressure P, 

= (p + P)u^u u -Pg^. (5) 
The equations describing the background evolution are 

ft + m{ P + p) = o, (6) 



o kS 2 u) ( ip' \ u>' 
dip b \ ip J <p 



tp" + 2Hip' 



kS 2 



2w + 3 



(P-3P), 



(7) 



(8) 



where the prime denotes derivative with respect to con- 
formal time T), S is the scale factor, and Tt = S'/S. Gen- 
eral relativity is recovered in the limits 



LU — > oo, tp 



0, P" 



0. 



(9) 



To recover the value Newton's gravitational constant 
today which is determined by Cavendish type experi- 
ments, we also require that the present day value of p 
is given by 



Po 



2w + 4 
2w + 3' 



(10) 



III. PERTURBATION THEORY 



A. The 1+3 covariant decomposition 



ip = G4>, 



(2) 



where G is the Newtonian gravitational constant mea- 
sured today. The Einstein field equations are then gen- 



The main idea of the 1 + 3 decomposition is to make 
space-time splits of physical quantities with respect to 
the 4-velocity u a of an observer. There are many possi- 
ble choices of the frame, for example, the CMB frame in 
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which the dipole of CMB anisotropy vanishes, or the lo- 
cal rest-frame of the matter. These frames are generally 
assumed to coincide when averaged on sufficiently large 
scale. Here it will be convenient to choose u a to coincide 
with the velocity of the CDM component, since u a is then 
geodesic, and acceleration vanishes. From the 4-velocity 
it a , we could construct a projection tensor h a b into the 
space perpendicular to u a (the instantaneous rest space 
of observers whose 4-velocity is u a ): 

hab = gab ~ U a Ub, (11) 

where g a b is the metric of the spacetime. Since h a b is 
a projection tensor, it can be used to obtain covariant 
tensor perpendicular to u a , and it satisfies 

h a b = h {ab) , u a h ab = 0, h c a h cb = h ab , h a a = 3. (12) 

With the timelike 4-velocity u a and its tensor counter- 
part hab, one can decompose a spacetime quantity into 
irreducible timelike and spacelike parts. For example, we 
can use u a to define the covariant time derivatives of a 
tensor T b - C d ,.. e : 

ft c e = u a V a Tt.t, (13) 

furthermore, we can exploit the projection tensor h ab to 
define a spatial covariant derivative D a which returns a 
tensor which is orthogonal to u a on every index: 

D a T b - c d ... e = h a p h b . . . h c s h\ . . . //,"V"7 ' Y..„. (14) 

If the velocity field u a has vanishing vorticity, D a reduces 
to the covariant derivative in the hypersurfaces orthog- 
onal to u a . The projected symmetric tracefree (PSTF) 
parts of vectors and rank-2 tensors are 

V<„> = K b V b , (15) 



T(ab) — h( a c hb) d T c d = h( a c h b ) d T cd — - h cd T cd h a b- (16) 

One can also define a volume element for the observer's 
instantaneous rest space: 

£abc = VabcdU d — £[ a bc\ j (17) 

where r\ a bcd is the 4-dimensional volume element (i] a bcd — 

V[abcd], ?7oi23 = -V\g~\)- Note that D c h a b = = D a £ bc d- 
The skew part of a projected rank-2 tensor is spatially 
dual to the projected vector T a — \e a bcT^ bc \ and any 
projected second-rank tensor has the irreducible covari- 
ant decomposition 

Tab = g Th a b + £abcT C + T( a b) , (18) 

where T = T cd h cd is the spatial trace. In the 1+3 covari- 
ant formalism, all quantities are either scalars, projected 



vectors or PSTF tensors. The covariant decomposition 
of velocity gradient are 

V a u b = D a u b + u a A bl (19) 



D a U b = Uab + <?ab + -^Ohab, (20) 

where a a b = D(a M &) is the shear tensor which satisfies 

Oab = <7(ab), <?a = and U a <7 a b = 0; U a b = D[ Q Ufc] 

is the vorticity tensor, which satisfies io a b — ^>[ab] and 
u a u! a b — 0. One can also define the vorticity vec- 
tor ui a — E ab c^ bc /2 (with Lu a b — Eabc^)- The scalar 
9 = V a u a = D a u a = 3H is the volume expansion rate, H 
is the local Hubble parameter; and A a = U WbU a = u a is 
the acceleration, which satisfies u a A a — 0. We note that 
the tensor D a u b describes the relative motion of neigh- 
bouring observers. The volume scalar 9 determines the 
average separation between two neighbouring observers. 
The effect of the vorticity is to change the orientation 
of a given fluid element without modifying its volume or 
shape, therefore it describes the rotation of matter flow. 
Finally, the shear describes the distortion of matter flow, 
it changes the shape while leaving the volume unaffected 

Gauge-invariant quantities can be constructed from 
scalar variables by taking their projected gradients. The 
comoving fractional projected gradient of the density 
field of a species i is the key quantity of covariant 
method 

X« = J^pW (21) 

which describes the density variation between two neigh- 
bouring fundamental observers. The comoving spatial 
gradient of the expansion rate orthogonal to the fluid 
flow is 

Z a = SD a 9, (22) 

which describes perturbations in the expansion. These 
quantities are in principle observable, characterizing in- 
homogeneity in a covariant way, and vanishes in the 
FLRW limit. 

The matter stress-energy tensor T^T 1 ' can be decom- 
posed irreducibly with respect to u a as follows: 

T ir ) = P u aUb + 2u (Q 9 b ) - Phab + n ab , (23) 

where p = T^u a u b is the density of matter measured by 

an observer moving with 4-velocity u a , q a = h b l T b " l ^u c is 
the relativistic momentum density or heat (i.e. energy) 
flux and is orthogonal to w", P e -h ab T^ b Tl) /3 is the 
isotropic pressure, and the projected symmetric traceless 
tensor n a b = ^<™6> is the anisotropic stress, which is 
also orthogonal to u a . The quantities p, P, q ai n a b are 
referred to as dynamical quantities and a a b, ^>ab, 0, A a 
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as kinematical quantities. In the FLRW limit, isotropy 
restricts T^™' to the perfect-fluid form, so the heat flux 
q a and anisotropic stress Tr a b must vanish. 

The remaining first-order gauge-invariant variables 
which we need are derived from the Weyl tensor C a bcd, 
which is associated to the long-range gravitational field 
and vanishes in an exact FLRW universe due to isotropy. 
In analogy to the electromagnetic field, the Weyl tensor 
can be split into electric and magnetic parts, denoted 
by E a t and B a b respectively. They are both symmetric 
traceless tensors and orthogonal to u a , 



E a b = C ac bd.u c u d — E <a b > , 
B ab = -*C acbd u c u d = -\e a ef C, 



efbdl 



— B <a 



<ab>- 



(24) 



(25) 



Here * denotes the dual, *C ac bd — ^Vac @efbd- 

For the radiation field, we can make a 1+3 covariant 
decomposition of the photon 4-momentum as 



p a = E(u a + e°), 



(26) 



where E = p a u a is the energy of the photon. e a describes 
the propagation direction of photon in the instantaneous 
rest space of the observer. The observer can introduce a 
pair of orthogonal polarization vectors (ei) a and (e2) a , 
which are perpendicular to u a and e°, to form a right- 
handed orthonormal tetrad {u a ,(ei) a ,(e2) a ,e a } at the 
observation point. The (screen) projection tensor is de- 
fined as 



"Hab = g a b - u a u b + e a e b , 



(27) 



which is perpendicular to both u a and e a , and satisfies 
H%(e } ) b = (ei)». _ 

Using the polarization basis vectors, the observer can 
decompose an arbitrary radiation field into Stokes pa- 
rameters I(E,e a ), Q(E,e a ), U(E,e a ) and V(E,e a )^. 
Therefore one can introduce a second-rank transverse po- 
larization tensor P a b(E, e c ) 



Pab(e l ) a (e J ) b 



If I + Q U + V 
2\U-V I-Q 



(28) 



for i and j = 1,2, and we have omitted the arguments 
E and e a . P a b oc E 3 TL a Hffcd, where fd is photon dis- 
tribution function. Decomposing P a b(E,e d ) into its irre- 
ducible components, one obtains 



Pab(E,e d ) 



I(E,e d )Hab + Vab(E,e d ) 



+ -V(E,e d )e a bce c 



(29) 



where the linear polarization tensor V a b(E, e d ) satisfies 
VaMYietf = ~ I g U Q Y (30) 



It is convenient to define energy-integrated multipole 
for total intensity brightness and the electric part of the 
linear polarization: 

I Al = J AE j dfi I{E,e c )e <Al> , (31) 
Ca, = M t 2 dE J dn e {Al _ 2 V ai _ iai )(E,e c ) (32) 



where eA, 



e a e b e c ...ei and M t 



^2I(J -!)/[(! +1)G + 2)]. 



B. The linearized perturbation equations 

In the 1+3 covariant approach, the fundamental quan- 
tities are not the metric, which is gauge-dependent, but 
the kinematic quantities of the fluid, namely the shear 
a a b, the vorticity uj a b, the volume expansion rate 9 and 
the acceleration A a , the energy- momentum of matter 
and the gravito-electromagnetic parts of the Weyl ten- 
sor. The fundamental equations governing these quan- 
tities are the Bianchi identities and the Ricci identities. 
The Riemann tensor in these equations is expressed in 
terms of E a b, B a b and the Ricci tensor R a b- The modified 
Einstein equation connects the Ricci tensor to the mat- 
ter energy-momentum tensor. In the following, we have 
linearized all the perturbation equations. We should also 
note that the definitions of covariant variables do not as- 
sume any linearization, and exact equations can be found 
for their evolution. 

The first set of equations are derived from the Ricci 
identities for the vector field u a , i.e. 



2V[ a V b ]U c = R ab cdU d . 



(33) 

Substituting the 4-velocity gradient (fl"9]) and the decom- 
position of the Riemann tensor, and separating out the 
time-like projected part into the trace, the symmetric 
trace-free and the skew symmetric parts, we obtain three 
propagation equations. The first propagation equation is 
the Raychaudhuri equation, 

e + \e 2 -D a u a + ( P + 3P) + 

I (2uj^ + -D a D a <p + 6^ + 3^) = , (34) 
z \ ip z w If tp / 

which is the key equation of gravitational collapse, ac- 
counting for the time evolution of 9. The second is the 
vorticity propagation equation, 

2 

Uab - D[ a U b ] + -9u a b = . (35) 

The last one is the shear propagation equation, 

i 2 a n • i z? i K7Tab 

°"<a6> + o t ' cr afc — U <a Ub> + Hiab + 77 

3 2 If 

+ -LD < bD a> ^+^(7ab = : (36) 

Zip lip 
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which describes the evolution of kinematical anisotropies. 
It shows that the tidal gravitational field E ab and the 
anisotropic stress ir a b would induce shear directly, and 
the shear will change the spatial inhomogeneity of the 
expansion through the constraint equations (|37p ■ 

The propagation equations are complemented by three 
constraint equations, which are spacelike components of 
Eg. ([55)) . The first is the shear constraint, 



D b u ab + D b cr ab - %D a d - —q a - LU^rD a ip 
- — (D a <p)- - —it a = 0, 

ip ip 



(37) 



which shows the relation between the momentum flux 
q a , the shear a a b and the spatial inhomogeneity of the 
expansion. The second constraint equation is the vortic- 
ity divergence identity, 



D c (e abc iu ab ) = 0. 



(38) 



and the second propagation is the £?-equation 



B a b + (0 + —)B ab — 

Zip 



D c E d(a + —D c ir d(a 



l-D c D d D {a ip - ±-D c D^p> h d(a 
= . 



(43) 



This pair of equations for electric and magnetic parts 
of the Weyl tensor would give rise wavelike behavior for 
its propagation: if we take the time derivative of the in- 
equation, commuting the time and spatial derivatives of 
B term and substituting from the B-equation to elimi- 
nate B, we would obtain a E term and a double spatial 
derivatives term, which together give the wave opera- 
tor acting on E; similarly we can obtain a wave equa- 
tion for B by taking time derivative of the i?-equation. 
These waves are also subjected to two constraint equa- 
tions, which emerge from the spacelike component of the 
decomposed Eq. (|4"T|) . The first constraint is 



The last one is the B ab equation, 



B, 



b + (D c uj d(a + D c a d(a )r] b)c d u e = , 



(39) 



which shows that the magnetic Weyl tensor can be con- 
structed from the vorticity tensor and the shear tensor. 
With this last equation B a b may be eliminated from some 
equations in favor of the vorticity and the shear. 

So far we have only discussed propagation and con- 
straint equations for the kinematic quantities. The sec- 
ond set of equations arises from the Bianchi identities of 
the Riemann tensor, 



V[ e i? C( j]afc — 



(40) 



which gives constraint on the curvature tensor and leads 
to the Bianchi identities for Weyl tensor after contracting 
once, 



V d C a6c d = V [h i? a]c + - g c[b V a] R. (41) 



D b E ab - ^(2D aP + 20q a + 3D b n ab ) + ^p^f- 
b<p 3 ip z 



K if 



1 ^ <P A L 



> ^ qa ~ ( 3 + 2V [ 3 J% + {DaV) ' + " = -( 44 ) 



This is the div E equation, with the source term given by 
the spatial gradient of energy density. It can be regarded 
as a vector analogue of the Newtonian Poisson equation, 
and shows that the scalar modes will result in a non-zero 
divergence of E ab , and hence a non-zero gravitational E- 
field. The second constraint equation is 



D b B, 



ab -7Ti(P + P^ab^Ucd + Vabc d U b D° q d ] 
Zip 

(^-^D,D^-^ + ^u b . cd + 
ip 3p 6 ip tp 



Va bcd u b L^D c D d ip + -(D c D d ^y 
\ ip z ip 



3ip 



ip 

D c D d ip 



= 



(45) 



The 1 + 3 splitting of the once contracted Bianchi iden- 
tities leads to two propagation and two constraint equa- 
tions which are similar in form to the Maxwell field equa- 
tions in an expanding universe, governing the evolution of 
the long range gravitational field. The first propagation 
equation is the ^-equation, 

E ab + 0E ab + D c B d(a r) b) d u e + ^-[3(p + P)a ab + 

3D <a qb> ~ 37T ab - O-Kab] + \<7 a b{u + - 

2 2 <pr 



1 Tab 



1 , 3.<f 



~^ Dfl D^ + ~(. + -)^D <a D b> ip 



1 ip 3 ip 

+ ~—E a b + -K^TTab = , 

2 ip 4 ip A 



(42) 



This is the div B equation, with the fluid vorticity serving 
as source term. It shows that the vector modes will re- 
sult in non-zero divergence of B a b, and hence a non-zero 
gravitational B-field. The above equations are remark- 
ably similar to the Maxwell equations of the electromag- 
netism, so we have chosen to use E a b and B a b as the 
symbols. 

The last set of equations arises from the twice- 
contracted Bianchi identities. Projecting parallel and or- 
thogonal to u a , we obtain two propagation equations, 



6{p + P) + D a q a = 



(46) 



q a + ^0q a + (p + P)u a + D b ir ab - D a P = , (47) 
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respectively. For perfect fluids, these reduce to 

p + 6(p + P)=0, (48) 



+ P)ii a - D a P = 



(49) 



which are the energy conservation equation and momen- 
tum conservation equation respectively. 

The background field equation for Brans-Dicke field 
is given in Eq.®. The first order covariant and gauge- 
invariant perturbation variable of the Brans-Dicke field is 
defined as the spatial derivative of the Brans-Dicke field, 



V a = SD a ip . 



(50) 



Taking the covariant spatial derivative of Eq.©, com- 
muting the spatial and time derivatives of V term, we 
could obtain the first order perturbation equation for 
Brans-Dicke field after linearization, 



V a ' + 2HV' a + SZ a cp' + S 2 D a D b V b 
kS 2 



3 + 2w ^ 



(51) 



The Eg. ((57)) is also the generalized Friedmann equation, 
showing how the matter tensor determines the 3-space 
average curvature. 

The last first-order covariant and gauge-invariant vari- 
ables can be obtain from the spatial derivative of the 
projected Ricci scalar, 



Va = ^SD a ^R , 

after a tedious calculation, we obtain 

pX a pV a 1 <p' 
Va = K k — 5 „(2rt H )Z a 



(58) 



•r 



if 2 s 



1 



1 p 1 ' 



o (p & p <p 



w ^ v 1 n n v 3 -w2 V a 

'T^ - T^a P>aP>vV ~ SlL 

b z ifr p bp 



(59) 



C. Mode expansion in spherical harmonics 



where the upper index (i) labels the particle species. 

In the absence of rotation, uj a b — 0, one can define 
a global 3-dimensional spacelike hypersurfaces that are 
everywhere orthogonal to u a . This 3-surfaces is meshed 
by the instantaneous rest space of comoving observers. 
The geometry of the hypersurfaces is determined by the 
3-Riemann tensor defined by 



[D a ,D b ]u c = ®R a , 



bdcU 



(52) 



which is similar to the definition of Riemann tensor R a bdc 
but with a conventional opposite sign. The relationship 
between (3) i? fcdc and R a bdc is 

{Z) Rabcd = -ha q h b S hJhd P R qs fp-VacVbd + V ad V bc 



= ^R 



ab] [cd] 



(53) 



where v ab = Df,w a is the relative flow tensor between two 
neighbouring observers. In analogy to 4-dimension, the 
projected Ricci tensor and Ricci scalar arc defined by 



and 



^Rab = ^Racbdh Cd = ^ R C acb , 



^R = ( ^R nh h ab . 



(54) 



lb a . (55) 
The ^R a b is determined by the Gauss-Codacci formula 

{3) Rab = \ {5) RK b - h<i ab - |tt qo + E ab , (56) 



where 



^R = 2( Kp - 1 -9 2 ) 



(57) 



In the linear perturbation theory it is convenient to ex- 
pand the 0(e) variables in harmonic modes, since it splits 
the perturbations into scalar, vector or tensor modes and 
decouples the temporal and spatial dependencies of the 
1+3 equations. This converts the constraint equations 
into algebraic relations and the propagation equations 
into ordinary differential equations along the flow lines. 
In this paper we focus on the scalar and tensor pertur- 
bation modes, since the vector modes would decay in 
an expanding universe in the absence of sources such as 
topological defects. 



1. Scalar mode 

For scalar perturbations we expand in the scalar eigen- 
functions Q'*' of the generalized Helmholtz equation 



S 2 D a D a Q 



(fe) 



k 2 Q 



(k) 



(60) 



at zero order. They are defined so as to be constant along 
flow lines, i.e. independent of proper time = O(e), 
and orthogonal to the fluid 4- velocity u a . 

For the I— th multipoles of the radiation anisotropy and 
polarization, we expand in the rank-^ PSTF tensor, Qj{ , 
derived from the scalar harmonics with 



D< ai ,..D ai >Q 



(fe) 



(61) 



where the index notation Ai denotes the index string 

(k) 

a\...ai. The recursion relation for the Q A , 
S 



: D <a l Q ( A ^_ 1> 



(62) 
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follows directly from Eq. (j61~j) . The factor of (S/k) 1 in the 
definition of the ensures that = Q at zero-order. 

(k) 

The Q A also satisfies some other zero-order properties, 



)(*) 



= 0, 



/'""<<>:':.„..„.,. (63) 



We also have the following differential relations which 
can be derived from Eqs. (|6"0)) and (|6"2")l : 



D 2 Qil 



I 



21-1 S 

k 2 ' 



I -{I 2 - 1) 



K 



S 2 



1-^ + 1)5 



)(fc) 

'ai...ai 



QS..a,(64) 

(65) 



Now we can expand the gauge-invariant variable in the 
following dimensionless harmonic coefficients: 



z n = 



q (i) 



E • 

fc ^ 

fc 

EW> 



\ - * n (J 



(fc) 
6 ' 



A, 



E 



E fc ^ fc) > 
E %^ 



Wn(fc) 

A, ' 



k 



(OnW 

A; ' 



(66) 

(67) 
(68) 
(69) 
(70) 

(71) 

(72) 

(73) 
(74) 

(75) 
(76) 
(77) 



where the upper index (i) labels the particle species. The 
scalar expansion coefficients, such as X„ , are first-order 
gauge-invariant variables, and their spatial gradients are 
second-order, for example D a X^ = 0(2). In the co- 
variant and gauge-invariant approach, we characterize 
scalar perturbations by requiring that the vorticity and 
the magnetic part of the Weyl tensor be at most second- 
order. Demanding u> a b = 0(2) ensures that density gra- 
dients are not from kinematic effects due to vorticity, and 



setting B a b = 0(2) ensures that gravitational waves are 
excluded to the first order. 

To obtain the scalar equations for the scalar expan- 
sion coefficients, one could substitute the harmonic ex- 
pansions of the covariant variables into the propagation 
and constraint equations given in the section above. Here 
we will consider only the adiabatic modes. For the (i) 
fluid, 



jjap{i) _ c (i)2 



(78) 



where Cg is the adiabatic sound speed of the (i) fluid. 
For the spatial gradients of the total density X k , we find 

k 



— Ep (i) 4 i) lc W2 -- 



E^ } J -3^(1 + ^)^=0. 



(1 + -)z k 
P 



(79) 



For the individual fluid of the (i) species, the propagation 
equation satisfies 



-3W(1 + — )W k = . 



■k[(l 



p(0, 
1M 



)Z k 



p 



For the heat fluxes, we have 



(80) 



(1 + -r^)kW k 



3 fc(1 k 2 



fccW 2 4° =0. 



(81) 



The heat flux for each fluid component is often given by 
q k — (p^ + pt'')i)[' , so we can derive the propagation 
equations for v k from Eq. (|81|) . 

We also can obtain the time evolution of the spatial 
gradient of the expansion 



z; ~^ { :m' m 2 fc 2 ) • ^E :l 



I kS 2 
~k~2Lp 



fc 2 n , V" u/ , if 12 

— ] + 4w^V' k + 3-^ + kZ k ^- + W k 4w-^V 

If Lp ip ip (p 

rJ> ,J ,J 



if) 

h 1 



H } ) P () x k +-\v k [-^—-i--—{ P + 

3P) 

+6— -3ft— ) +3— Wl) = . (82) 

ip LP LP J 
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Substituting the covariant harmonic expansion into 
Eq. (|59p , and then taking the time derivative of this equa- 
tion, we obtain the evolution of the spatial gradient of the 
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3-curvature scalar: 



pressed as 



k 2 r,' k = -X k [S 2 n(p + 3P)- + S 2 k P ^} 

p ip z 



S 2 n{p + 3P)^ + ls 2 K (p - P)4 
p z 2 p 3 



in" m' 2 



m' 



3W 
2 ] ip 3 

2 V 



ip' 3 „ w' 
<p A 

tp' 2 p' 
p- 3 p z 



W fc (2w^--2w^ 



3H— + 3W^-) + W' k (uj 



12 



p 2 



m—) 



12 i 

kZ k (2H' + £--?—)- kZ' k (2H + t.) 



(83) 



P X k = pWjSf^ + pMlf 1 + P ( b 'xf + p (c) 4 c) (8J 



P9fc 



+P (C) 4 C) ■ 



(6) 
A- 



(89) 
(90) 



2. Tensor mode 

For tensor modes, we expand the first order pertur- 
bation variables in the rank-2, zero-order PSTF tensor 
eigenfunctions Q^j) G f the comoving Laplacian, 

S 2 D C D C QW =k 2 Q^. 



(91) 



As mentioned by Ref . [511 ] , solving the propagation equa- 
tion of r] k avoids the numerical instability problem in 
isocurvature modes when we work in the CDM frame. 

From the shear propagation equation (|36p . the propa- 
gation equation for Ok becomes 



a' k + Ha k - kW k + k<P k + — — pn k + k— 



k 2p' 



J 2cp 



It 

2 p 



<Tk =0 



(84) 



Similar to the case of scalar modes, this equation holds 
at the zero-order. The tensor harmonics are transverse, 
orthogonal to u a , and constant along the integral curves 
of u a : 



o, 0$ = 



(92) 



They can also be classified as having electric parity (de- 
noted by Q^}) or magnetic parity (denoted by Q^)- 
These two parity harmonics are related by a curl: 



From the Div E equation (|44|) . we could obtain the <fr k 
equation, 



k 3 

37i up 

S p 



^[X fe + (1-3^M- 



k 

k p' 



3 1 p' 
2SP 2 







(85) 



The algebraic equation of a k can be derived from the 
shear constraint equation (|3"7) 



<7fe(l 
1 



4»i 



211 

k p 



■ V fe 



(V£ - WV fe ) + ^VF fe = 



(86) 



curl 



(fc) 



3if 



(93) 
(94) 



For tensor mode perturbations, the vorticity and all 
gauge-invariant vectors vanish at the first order, i.e. 
ujab, X a ,Z a ,q a ,Aa,Va, r\ a all equal to zero [H, The 
rest rank-2 gauge-invariant tensors are constrained to be 
transverse: 

{3) V a E ab = 0, {3) V a B ab = 0, 

< 3 ) W„6 = 0, (3) W ah = 0. (95) 

And they can be expanded in electric and magnetic parity 
tensor harmonics: 



From the first-order perturbation equation for the Brans- 
Dicke field (p)Tj) , we could derive the quadratic differential 
equation of Vk 



V'k' 



f 2HV' k + kZ k p>' + k 2 V k 
kS 2 



3 + 2w 



2(l-3cW 2 )pWjfW 



(87) 



The variables X k , qk and ilk (without upper index (i) 
) refer to variables of the total matter, and can be ex- 



E„ 



Bab = 



Cab = 



E 

k 

E 

k 

* ^ k 



-<p( E kQab + E kQ K ab) 



^{B k Q a k b ] + B k Q a ' 



(a k Q a y +v k Q a l') 



(96) 
(97) 
(98) 



t2 = ^EfrW+^W). 09) 
fc 
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Substituting these into equations in section IIII B[ we 
obtain the E k and <r k propagation equations: 

k 2 E' k + k 2 E k (H + - k 3 (l + 3^)a k + 

— - (p + P fcdfe + — ^ (uj + -)ka k + H + 

zip ztp z 2 tp 

3nS 2 Pir k 3^ 2 1 2 , 

W + -:^kS pir k - —kS pir k = ,(100) 

2 (p 4ip z 2ip 

a k = -Ha k - kE k - —— -— . (101) 

2k tp 2 ip 

IV. NUMERICAL IMPLEMENTATION 

We carry out our numerical study by modifying the 
CAMB code. The original CAMB code, written by 
Antony Lewis and Anthony Challinor[53j, is a For- 
tan 90 program which calculates CMB anisotropies in 
the standard Einstein general relativity, by solving the 
Boltzmann-Einstein equations for various components in 
the Universe. Most of the equations to be solved are 
in the file equations. f90, which can be modified conve- 
niently. The background evolution equation dr/da is 
written in function dtauda, and it can be modified for 
different background. The Boltzmann-Einstein equation 
group is listed in the functions fderivs (scalar mode for 
flat Universe), fderivst (tense mode for flat Universe), 
derivs (scalar mode for non-flat Universe) and derivst 
(tense mode for non-flat Universe). This equation group 
includes the propagation equations of scalar factor S, 
the 3-Ricci scalar perturbation 77, the cold dark mat- 
ter perturbation X c , the baryon perturbations Xf, and 
Vb, photon multipole moments, and neutrino multipole 
moments in the covariant approach. The CAMB code 
uses the Runge-Kutta method (subroutine dverk in file 
subroutines. f90) to solve these equations. To speed up 
the calculation, the line-of-sight integration method first 
developed by Seljak and Zaldarriaga[54| is used: the dif- 
ferential equation for photon temperature perturbation 
is integrated along the l.o.s. to obtain ST/T. The mul- 
tipoles today is a definite integral of source term multi- 
plied by the spherical Bessel functions from early time 
to today. The source term of scalar perturbation at a 
given time for a given wavenumber is encoded in sub- 
routine output. The subroutine evolves the perturbation 
equations and does the integration in cmbmain.f90 file. 
The main routine for running CAMB is wrapped in file 
camb.f90. Using these equations, we modify the code 
for calculation in the Brans-Dicke theory. The most im- 
portant three parts of modifications are: the background 
evolution, the Boltzmann-Einstein differential equations, 
and the source term in the line-of-sight integration. 

For the background evolution, we implement the pro- 
cedure described in the appendix of Ref. [15]. To satisfy 
the end point condition Eq. (|10p , we start from an epoch 



which is deemed early enough. We then evolve the model 
forwards (to avoid numerical instability we do not evolve 
backwards) to obtain the tp value today. The procedure 
is repeated with a Brent algorithm (see e.g. [55|) to find 
the initial value of tp at that epoch. In doing this we 
set <p' = and V k = V' k = at the initial point. The 
initial condition ip' — can be justified by Eq.©: in the 
radiation dominated era, the R.H.S. of Eq.©, p — 3P is 
negligible compared with other terms, then 

tp' = ci +c 2 S- 2 . (102) 

This mean that any initial velocity quickly dies out in 
a few Hubble times and approaches a terminal velocity 
Ci, this velocity is constrained by nucleosynthesis, so it 
should be very small. The initial condition V k = V' k = 
is the simplest choice which matches the requirement of 
Eq. (|87p . Initial perturbations in tp are damped during 
the radiation dominated era, so the choice of the initial 
condition of V k have little impact on CMB anisotropy in 
the adiabatic perturbation case. 

To realize the background evolution described above, 
we write a separate module. The function of this module 
is that, for a given value of tp today which is determined 
by Brans-Dicke parameter w, we first find out the initial 
value of ip at sufficiently early time which can evolve the 
given value of tp today, and then we could calculate tp 
and tp' at each scale factor S and store them into arrays 
for interpolation in subsequent process. Therefore, if one 
want to use tp and tp' in the code, just simply use this 
module first. 

To be consistent with modified Friedmann equation 
([7]), in the code we define the critical density as 

Per = ^H 2 , (103) 

K 

where Hq is the hubble parameter today and tp is given 
in Eq. flTJ]) . This definition differs from the conventional 
one by an additional factor tp. The definition of the 
fractional density is the same as the traditional one: 
= Pq I Per ■ Because tp' approximately vanishes to- 
day (c.f. Fig.©), from Eq.© we find n to tai ^ 1 for 
the flat geometry. This definition is convenient in study- 
ing the non-flat universe. We also should note that the 
difference with the traditional one is very small, in most 
case, less than 1%, because tp = 1.001 when w — 50. 

In this work we adopt the cosmological constant as 
dark energy, this is equivalent to set the potential of the 
Brans-Dicke field to a constant. The more general case of 

extended quintessence [5H Unl IHS M, m, El win be 

dealt with in future studies. Below, we adopt the ACDM 
model with Einstein gravity which best fit the WMAP 
five-year data[63| as our fiducial model, i.e. £Ia = 0.721, 
n b = 0.0462, fi c = 0.233, H = 70.1 km s" 1 Mpc" 1 , 
n s = 0.960, = 2.457 x 10" 9 at k = 0.002 Mpc" 1 , 
Zrcion = 10.8 and the equation state of dark energy w — 
-1. 

For the the Boltzmann-Einstein differential equations, 
we modified the scale factor evolution equation and r\ k 
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propagation equation according to the Eqs. |J7J| and (|83|) 
respectively in functions fderivs and fderivst in equa- 
tions. f90. Some other complementary equations, such as 
the constraint equations, have also been modified corre- 
spondingly. 

To speed up the calculation, the CAMB code integrates 
the system of differential equations by using the line-of- 
sight integration method, first developed by Seljak and 
Zaldarriaga for the CMBFAST code [II]. In this method, 
the multipole moment of photon intensity 1^ could be 
express as [42|, [64| , 



= 4 



k 



dr) Se T \ [ — cjfc + -n e cr T K 2 ( 



«=-100 
w=400 
«=-400 ■ 



^1.00- 



10" 4 

s 



+<7 T v k ——^{x) 
kr ax 



Ik 1 (o) 

3 S 4 e 1 k 



3 $ ' W+ fc¥ dx 2 



(104) 



where 770 is the conformal time today, x = (ry — 

r = 1/ and r is the zero-order optical depth back 

to x. Here, 



(Z-i/)! a; 17 : 



(105) 



(1 



are the ultra-spherical Bessel functions, k 2 
3K /k 2 ) 1 / 2 , and £^ is the quadrupole of the E-like polar- 
ization of the CMB photons. After integration by parts, 
one could eliminate the derivatives of ultra-spherical 
Bessel functions and write temperature anisotropies as a 
time integral over a geometrical term $^ (ir) and a source 
term: 



r(0 _ 
1 k — 



dr] $?(x) x S 



(106) 



where the source term is given by 
1 



S = 



12k 2 



K'2 



12ka k e~ T n 2 + 24ka' k g(r))K 2 



12 ka k g'(ii)n 2 + 3g"( V )( k + 6g'(v)C k + 3g(v)( k 
+ 12 kK 2 g'(r))v k + 12 kK 2 g{rj)v' k + 4k 3 a k e~ T K 2 H 



k 2 g{v)(k -4k 3 e- T Z kK2 + 3k 2 g( V )li 0) K2 



in which 



Km 



9,(2) 

= n e <7TSe~ 



, (107) 



(108) 
(109) 



g(rj) is the visibility function. Using the first order deriva- 
tive perturbation equations described in Section IIIIl <r' k 



FIG. 1: The time evolution of the Brans-Dicke field tp. 



V 




FIG. 2: The evolution of time derivative of the Brans-Dicke 



and Cfc m the source terms could be further expanded 
to the zeroth and first order derivative terms which are 
expressed in variables used in the output subroutine in 
the CAMB code. 



V. RESULTS 

In Fig[V] we show the time evolution of the Brans- 
Dicke field tp. For models with u> > 0, the value of tp 
increases with time, whereas for models with u> < 0, tp 
decreases with time. During the radiation dominated 
era, the variation of ip is very small, almost zero. When 
entering the matter dominated epoch, tp begins to in- 
crease or decrease. After dark energy begin to dominate, 
tp changes more rapidly. The time evolution of tp' is plot- 
ted in Fig[V] I tp' I reaches a terminal velocity in the ra- 
diation dominated era, and then begin to decay in the 
matter dominated epoch. In the dark energy dominated 
age, it will increase again, and its present day value for 
this model is of the order 10~ 6 . 
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FIG. 3: The time evolution of the effective Newtonian gravi- 
tational coupling G e ff. 



The effective Newtonian gravitational coupling G e / / is 
the inverse of tp in the unit of G. The time evolutions of 
G e ff are shown in Fig|V] We can see that G e // changes 
rapidly at low redshift, so it may not be reliable to use the 
Type la supernovae (SNe la) data to constrain the Brans- 
Dicke theory: the Chandrasekhar mass Mch ^ G~ 3 / 2 , 
so the variation of the gravitational coupling G means 
that the peak luminosity of SNe, which is approxima- 
tively proportional to the Chandrasekhar mass, may also 
change, making it not reliable as a standard candle. 

The CMB temperature and polarization angular power 
spectra for the Brans-Dicke theories with u> = oo(i.e. 
general relativity) and ±100 are plotted in Fig|V] As 
can be seen, compared with the general relativity theory 
with the same cosmological parameters, both the loca- 
tion and height of the CMB acoustic peaks are changed. 
The Brans-Dicke model with a positive ui has broader 
and lower acoustic peaks for this set of parameters. As 
\oj\ increases, the difference in CMB angular spectra be- 
tween Brans-Dicke theory and general relativity dimin- 
ishes. The difference is more apparent at large I (small 
angular scale), so high resolution CMB data would be 
very useful in distinguishing the different models. From 
FigEl it is also very clear that the polarization spec- 
tra have a strong discriminating power. With the higher 
angular resolution and polarization data which we ex- 
pect in the nearby future, we should be able to lift the 
the degeneracy of parameters and place a more stringent 
constraint on the Brans-Dicke models. 

We compared the CMB angular power spectra results 
with those obtained with our CMBFAST code in the syn- 
chronous gauge [HI . We find that the difference in CMB 
power spectra is typically less than 1 percent [6(|, and is 
due primarily to the difference in the original (Einstein 
gravity) codes-for really making highly precise constraint 
on cosmological parameters with the CMB data, the pre- 
cision of the CMB Boltzmann needs to be further im- 
proved. The new code of course has better program ar- 
chitecture and runs faster. Particularly, if one calculate 




500 1000 1500 



FIG. 4: CMB temperature and polarization power spectra for 
Brans-Dicke theories with lu = oo, ±100 in scalar mode. 




FIG. 5: AC; = G{lo 
EE correlations. 



oo) -Ci(lj = 100) for TT, TE and 



dCi/dui, which reflects the impact of the gravity model on 
the CMB angular power spectrum, the results of the two 
code agree with each other at high precision, as shown in 
FigH The result on AG; = Ci{u = oo) - Ci{uo = 100) 
for TT, TE and EE correlations are very consistent in 
two codes, and the two curves are almost indiscernible in 
FigH 

We also plot the CMB temperature and polarization 
spectra yielded by tensor modes in Fig [6] The tensor- 
to-scalar ratio is set to 0.1. The primordial gravitational 
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FIG. 6: CMB temperature and polarization power spectra 
for the Brans-Dicke theories in tensor mode. The solid, dotted 
and dashed curves represent the Brans-Dicke model with ui = 
oo, 100 and —100 respectively. The tensor-to-scalar ratio R is 
set to 0.1 

wave produces large temperature fluctuations at the large 
scales, as well as a unique B mode polarization. In con- 
trast to scalar modes, when compared with the result 
of general relativity, the height of the peaks are higher 
for positive uj. Similar to the scalar mode, positive to 
shifts the peaks to smaller scale. At both the very large 
scales and very small scales, the differences in spectra be- 
tween the Brans-Dicke theory and the general relativity 
are very small, almost invisible, and yhe differences are 
only sensitive at I ~ 80. 

Fig[7]shows the impact of Brans-Dicke field on the mat- 
ter power spectra at z = 0. For u> — 100, the bend of the 
matter power spectrum occurs at short wavelengths, and 
there is thus more small-scale power, in agreement with 
the prediction of Ref. [65[ . 

VI. CONCLUSION AND SUMMARY 

Compared with Einstein's general relativity, there is an 
additional scalar field coupled with the Ricci scalar in the 
Brans-Dicke gravity, which makes the perturbation the- 
ory more complicated. With a covariant 1+3 approach, 
we have developed a full set of covariant and gauge- 
invariant formalism for calculating the cosmic microwave 
background temperature and polarization anisotropies in 
the Brans-Dicke gravity. Instead of using the components 
of metric as basic variables, the covariant formalism per- 
forms a 1+3 split of the Bianchi and Ricci identities, 
using the kinematic quantities, energy-momentum ten- 
sors of the fluid(s) and the gravito- electromagnetic parts 



10'F 




0.0001 0.0010 0.0100 0.1000 1.0000 

k [h Mpc~'] 

FIG. 7: The matter power spectra at z=0. The solid, dotted 
and dashed curves represent the Brans-Dicke model with u = 
oo, 100 and —100 respectively. 



of the Weyl tensor to study how perturbations evolve. 
Adopting covariantly defined, gauge-invariant variables 
throughout ensures that in our discussion the gauge am- 
biguities is avoided, and all variables had a clear, phys- 
ical interpretation. Since the definition of the covariant 
variables does not assume any linearization, exact equa- 
tions can be found for their evolution, which can then 
be linearized around the chosen background model. Fur- 
thermore, unified treatment of scalar, vector and tensor 
modes do not require decomposing the different modes 
from beginning as done in the metric method. A price 
we have to pay is that with this method the calculation 
is more complicated. 

We then calculate the CMB temperature and polar- 
ization spectra for the Brans-Dicke models using a mod- 
ified CAMB code. In this paper we consider both the 
scalar modes and the tensor modes in adiabatic initial 
condition, and adopt (fo — (2u> + 4)/(2u> + 3) at the cur- 
rent epoch and tp' = at early time as initial condition 
of the Brans-Dicke field. Compared with the general- 
relativistic model with the same cosmological parame- 
ters, both the amplitude and the width of the acoustic 
peaks are different in the Brans-Dicke models. We find 
that the small scale spectra and the polarization spectra 
will provide a sensitive and vigorous constraint on the 
different Brans-Dicke models in scalar mode. For tensor 
modes, the largest difference in CMB spectra for various 
Brans-Dicke models are located at I ~ 80. The struc- 
ture formation process in the Brans-Dicke theory is also 
studied. The matter power spectrum is shown in FigJT] 
For positive u case, the bend of the matter power spec- 
tra occurs at shorter wavelengths, and there is thus more 
small-scale power compared with the General Relativity 
case. 

Our numerical results confirmed the results obtained 
with particular guages (e.g. the synchronous gauge [15]). 
We also obtained for the first time the temperature and 
polarization spectra for tensor mode perturbations in the 
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Brans-Dicke theory. In the present code the speed and 
program structure are greatly improved, thus providing a 
more powerful and convenient tool for further studies. In 
paper II, we use the code and MCMC algorithm to derive 
the constraint on the Brans-Dicke parameter to with the 
latest CMB and LSS observational data. 

As the final remark, the covariant approach and cor- 
responding CMB code for the Brans-Dicke theory devel- 
oped in this paper, together with the synchronous gauge 
approach and corresponding code developed in the pre- 
vious paper [l5j, provide consistent, systematic and com- 
plete methods to study the Brans-Dicke theory. These 
methods and codes could be generalized to study more 
general scalar-tensor theory, as well as more complex ini- 
tial condition, we plan to carry out such generalization 
in subsequent studies. 
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